using System;
using System.Reflection;
using NBody;
using NBody.Analysis;
using OptionParser;

[assembly: AssemblyTitle ("NBodyInfo")]
[assembly: AssemblyVersion ("1.0.*")]
[assembly: AssemblyDescription ("Print information about the N-body system.")]
[assembly: AssemblyCopyright ("2005 Joseph D. MacMillan")]

public class InfoOptions : InputGetOptions
{
    [Option ("-f", "File for potential data")]
    public string PotFile;
    
    public InfoOptions()
    {
        PotFile = "";
    }
}

public class NBodyInfo
{
    public static void Main(string[] args)
    {
        InfoOptions opts = new InfoOptions();
        opts.ProcessArgs(args);
        
        NBodySystem s = NBodySystem.Read(opts.InFile);
        
        Console.Error.WriteLine("System:");
        Console.Error.WriteLine("  Number of particles = {0}", s.NumParts);
        Console.Error.WriteLine("  Time                = {0:f6}", s.Time);
        
        double m = 0.0;
        foreach (Particle p in s)
            m += p.Mass;
        
        Console.Error.WriteLine("  Mass of particles   = {0:f6}", s[0].Mass);
        Console.Error.WriteLine("  Total system mass   = {0:f6}", m);
        Console.Error.WriteLine("  Maximum radius      = {0:f6}", s.Rmax);
        Console.Error.WriteLine("  Minimum radius      = {0:f6}", s.Rmin);
        
        double xmin = Double.PositiveInfinity;
        double xmax = Double.NegativeInfinity;
        foreach (Particle p in s)
        {
            if (p.X > xmax) xmax = p.X;
            if (p.X < xmin) xmin = p.X;
        }
        
        Console.Error.WriteLine("  Maximum X position  = {0:f6}", xmax);
        Console.Error.WriteLine("  Minimum X position  = {0:f6}", xmin);
        
        double lambda_prime;
        double jx = 0, jy = 0, jz = 0;
        foreach (Particle p in s)
        {
            jx += p.Mass * p.Jx;
            jy += p.Mass * p.Jy;
            jz += p.Mass * p.Jz;
        }
        double jtot = Math.Sqrt(jx * jx + jy * jy + jz * jz);
        lambda_prime = jtot / Math.Sqrt(2.0 * m * m * m * s.Rmax);
        
        Console.Error.WriteLine("  Total angular momentum = {0:f6}", jtot);
        Console.Error.WriteLine("  Total specific angular momentum = {0:f6}", jtot / m);
        Console.Error.WriteLine("  Alt. Spin parameter    = {0:f6}", lambda_prime);
        
        if (opts.PotFile != "") // can calculate total energy
        {
            double Etot = 0;
            double lambda;
            
            Energy.SetPotentialFromData(ref s, DataFile.Table.Read(opts.PotFile));
            foreach (EnergyParticle p in s)
                Etot += p.Mass * (p.Potential + 0.5 * p.AbsoluteVelocity * p.AbsoluteVelocity);
            
            lambda = jtot * Math.Sqrt(Math.Abs(Etot)) / Math.Pow(m, 2.5);
            
            Console.Error.WriteLine("  Total Energy    = {0:f6}", Etot);
            Console.Error.WriteLine("  Spin  parameter = {0:f6}", lambda);
        }
        
        NBody.Cosmology.Model scdm = NBody.Cosmology.Model.SCDM;
        double z = scdm.Redshift(s.Time);
        Console.Error.WriteLine("  If cosmology is SCDM, then:");
        Console.Error.WriteLine("    z = {0}", z);
        Console.Error.WriteLine("    rho_b = {0}", scdm.Density(z));
        
        double v_avg = 0.0;
        foreach (Particle p in s)
        {
            v_avg += p.AbsoluteVelocity;
        }
        v_avg /= s.NumParts;
        Console.Error.WriteLine("Average velocity: {0:f6}", v_avg);
        Console.Error.WriteLine("Would take {0} time to travel across system", (xmax - xmin) / v_avg);
    }
}
